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The strong-shock, point-source solution and spherical isothermal distributions were used as initial condi- 
tions for a numerical integration of the differential equations of gas motion in Lagrangean form. The von 
Neumann-Richtmyer artificial viscosity was employed to avoid shock discontinuities. The solutions were 
carried from two thousand atmospheres to less than one-tenth atmospheres peak overpressure. Results 
include overpressure, density, particle velocity, and position as functions of time and space. The dynamic 
pressure, the positive and negative impulses of both dynamic pressure and static overpressure, positive and 
negative durations of pressure and velocity, and shock values of all quantities are also described for various 
times and radial distances. Analytical approximations to the numerical results are provided. 


I. INTRODUCTION 


HE problem of a spherical blast in air has been 
solved analytically for strong shocks in an ideal 
gas by J. von Neumann, and in the weak shock 
approximation by H. A. Bethe.! For intermediate 
regions it has been found necessary to resort to numeri- 
cal methods. The availability of high-speed computing 
machines has made possible the solution of this problem 
in the intermediate range with considerable accuracy. 
Two other groups are currently engaged in performing 
this work by independent methods. Under J. von 
Neumann and H. Goldstine at Princeton a shock 
fitting method has been employed. The differencing 
method of Peter Lax? is being used at New York 
University by S. Lowell. 

The method used in this paper is due to von Neumann 
and Richtmyer® and employs an artificial viscosity as 
a mechanism for avoiding shock-front discontinuities. 
Previously, T. S. Walton reported some results‘ using 
this method and an initial isobaric sphere of about 13 
atmospheres. 

The integrating process consists of the step-wise 
solution of difference equations which approximate 
the differential equations of motion of the gas. In 
order for such a procedure to be workable, however, a 
number of practical conditions must be satisfied. The 
differencing scheme must be stable, must offer reason- 
ably detailed results, must conserve numerical signifi- 
cance, and when put in the form of coded instructions 
for a high-speed computer, must be fast enough to 
reach desired solutions with a reasonable expenditure 
of machine time. - 


Il. METHOD FOR NUMERICAL INTEGRATION 
Equations of Motion 


The Lagrangean equations of motion are reduced to 
dimensionless parameters wherein pressure (P), density 


1 Shock Hydrodynamics and Blast Waves (Los Alamos Scientific 
Laboratory, 1944), Report AECD~2860. 

2P, D, Lax, Commun. Pure Appl. Math. (Institute of Mathe- 
matical Sciences, New York University) VII, 159-193 (1954). 

sj. ha Neumann and R. D. Richtmyer, J. Appl. Phys. 21, 232 
(1950). . 

47. S. Walton, Phys. Rev. 87, 910(A) (1952). 


(p), and velocity (#) are measured in units of ambient 
pressure (Po), density (po), and sound velocity (Co), 
respectively. In the following discussion, the pressure 
will frequently be expressed in atmospheres where 
one atmosphere is defined as equal to the pre-shock 
ambient pressure (Po). Where the expression over- 
pressure (AP) is used the reference is to the pressure in 
atmospheres in excess of the ambient pressure (Po), 
i.e., AP=P—1. We shall continue to speak of a velocity 
(u), however, when more properly we might refer to a 
mach number (u). The radial distance 7(ro,t), is expressed 
in energy-reduced dimensionless units (ro being 
Lagrangean distance, and / the time), such that 


AN=r/e and A=710/e, (1) 


where ¢ is a length expressing the energy and ambient 
pressure scaling: 


Etor 4 r® uv? 4nR? 
é=—=— (Bat) rdr————. (2) 
Po Povo 2/ ~  3(y—-1) 


Exot is the total blast energy and Fint is the specific 
internal energy. The subtracted term represents the 
pre-shock internal energy of the gas engulfed by the 
shock, and R is the shock radius. Time (¢) is defined in 
dimensionless units (7) such that 


r= tIto/e. 


The artificial viscosity (g), which acts like a pressure, 
is in units of the pre-shock ambient pressure (Po). 

In these units the Lagrangean equations of motion are 
written as follows: 


a 1 dp 2u du/dx 
=a or =n o(— +) (mass), (3) 
ax pr? ar dr» dr/dx 
ou —N A 
—=———(P+q) (momentum), (4) 
Or yy OX 
OP 14 
—=-—[yP+(y—-1)q] (energy), (5) 
Or par 
On 
“=, (6) 
Or 
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In these equations the Lagrangean variable is taken 
to be x=4(ro/6)?. 

In Eq. 5 (energy conservation) the internal energy for 
an ideal gas has been used. 


‘ie P po (7) 
p(v—1) Po 


Form of Artificial Viscosity 


An appropriate viscosity for the case of an outward 
moving spherical shock wave is the following: 
)® 


ov(v¥thsMy? du fdu |d 
4 3r 0x7 \ dx 


at 
Ox 

where Ax is the grid size and M is the number of grid 
zones in the shock front. The form of g chosen here is 
basically similar to that of von Neumann and 
Richtmyer, and satisfies similar conditions for the 
spherical case. However, there is no convenient steady- 
state solution by which to demonstrate the connection 
of variables across the shock, and the form chosen here 
is only asymptotically the same as a verified form for a 
plane problem. 

The fact that the plane form approximation is good 
for shock thicknesses which are small relative to the 
shock radius is borne out by the close agreement with a 
similarity solution (with g) investigated by R. Latter® 
in which he solves numerically the resulting ordinary 
differential equations for cases of interest to this 
problem. 

The particular form chosen for the viscosity has the 
advantage that it contributes nothing in regions of 
expansion, and is nonzero only in the compression 
phase of a shock. In Lagrangean coordinates this has 
the advantage of eliminating a spurious contribution 
near the origin where the positive velocity gradient is 
large. 


Initial Conditions 


Two general types of initial conditions were taken, 
(1) a point source, and (2) an isothermal sphere. A 
point-source case was run in the greatest length using 
the von Neumann strong-shock solution beginning at 
1600 atmos shock overpressure, and running down to 
less than 0.06 atmos. The point-source solution, 
starting at 199 atmos peak overpressure, was used to 
start the one problem carried out on the ORDVAC 
machine at the Ballistics Research Laboratories, 
Aberdeen, Maryland. The latter problem was run to 
nearly 0.1 atmos shock overpressure. The point-source 
solution was also used for starts at 473 and 818 atmos 
shock overpressure, and these problems were run down 
to less than 100 atmos peak overpressure. 

Three isothermal-sphere type problems have been 
run so far. Two of these begin with a hot isothermal 


5R, Latter, J. Appl. Phys. (to be published 1955). 
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sphere for which the density is the same inside and out. 
One began at 2002 atmos, the other at 121 atmos, 
overpressure. 

One cold isothermal sphere was also used, beginning 
with 121 atmos overpressure and a temperature inside 
equal to the outside temperature. 

These three problems were carried down to shock 
overpressures less than unity. 


Stability and Problem Checks 


Short plane-geometry problems were found useful 
in checking the validity of various computational 
tricks, and in studying the effect of certain changes in 
the difference equations. The effect of the size of the 
viscosity constant was also more easily studied in the 
plane case. Time increments smaller than that pre- 
scribed by the Courant condition (At= (pp)-*Ar) by a 
factor of 2 to 4 were found necessary for smooth results. 
In an attempt to reduce the required computing time 
the stable differencing scheme of Du Fort and Frankel® 
for diffusion-type equations was employed in place of 
explicitly carrying a viscosity quantity (g). Unfortun- 
ately such a scheme has some practical disadvantages. 
It requires carrying throughout a machine calculation 
sets of data for all space points for two different times. 
Furthermore, computing, changing time increments, 
and combining space points all become more tedious. 
Besides these disadvantages, additional terms must be 
introduced to correct for the excess energy introduced 
by the differencing scheme. On the other hand, the 
very general nature of the viscosity method, the ease 
of its applicability, and the precision with which it 
reproduces the Hugoniot conditions across a shock 
would seem to offset the more stringent time require- 
ment. Use of this method for nonideal gases is not 
considered in this paper, however. 

It is frequently convenient to use unequal zone sizes. 
For instance, the use of small zones through the shock 
front provides a sharp shock at very little cost in 
computing time. The use of such unequal zones was 
empirically validated in this problem by repeating 
calculations with quite different zone choices. 

The size of the time increments were automatically 
doubled whenever the stability conditions would allow 
it. Two conditions exist, one being the usual Courant 
condition and the other a diffusion-type condition 
imposed by the artificial viscosity 


Arg Ax/ Ne (Pp) *max 
(9) 


y 1 |d 
Ar<—( as) — 
4 Ox 


rq an 


: 
min 


The total energy in the blast wave must be conserved, 
but a check on the total energy is not a very sensitive 
test of the correctness of the results. Several machine 


6E. C. Du Fort and S. P. Frankel, Mathematical Tables and 
Other Aids to Computation, Vol. VII, p. 135 (1953). 
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errors were detected (and eliminated) which did not 
substantially effect the total energy. Here the system 
employed was to back off and re-run whenever sus- 
picious aberrations occurred and to re-run and com- 
pare the entire problem with a different zone spacing, 
different viscosity, different time increments, as well 
as some differences in the difference equations. 


Difference Equations 


The differential equations are approximated by the 
following difference equations: 


Ar(Ar")? 
uyrtt = uy? ——_—_—_—- 
(Ax) ty 
[Pav — Pa +gua* oy"), (10) 
AMM =D + "tH 1Az, (11) 


1—W 
pith=puin( ), (12) 
1+W 
where 
2(uy"44+-_1"*}) 
=Ar ( Se ee nee rere 
AA AL at a” 
unt — yy rtd 
APH —Ari*tH— ALA" 
y(y+1) s M3 
gyrtis =~) pry? Lepattt— art Pp, 
2 39 


for Uy_artt > uy"*3, 


Fic. 1. Radial dependence of the peak or shock overpressure. 
The solid curve represents the point-source solution. The dashed 
curves represent the results from two different initial isothermal 
spheres, one at 2000 atmos and the other at 121 atmos, both 
spheres were hot with equal density inside and out. The dotted 
curve represents the resulting shock overpressure from an initially 
cold isobaric sphere with temperature equal inside and out. The 
pressures are in atmos and the distances in units of (Eto¢/Po)t. 
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Fic. 2. Particle velocity or Mach number at the shock as a 
function of shock radius for the point-source solution. Velocity 
is in units of the pre-shock sonic velocity and the radius is in units 
of (Etot/Po)t. 


gattt=0 for uy tts ur}, 

y+1 

ig aia Pi” 

y— 

+2 (p44! — pry") q_y""? 
Piyt= » (15) 

pee 
——p1-3"— py" 
y-1 ° 


II. RESULTS OF NUMERICAL SOLUTIONS 
Peak Overpressure (P,) versus Shock Radius (2,) 


The strong-shock point-source solution provides that 
the shock overpressure (AP.) should depend on the 
inverse cube of the shock radius (A,). The problems 
that were begun with strong-shock, point-source values 
continued to obey the inverse cube law down to 10 
atmos at which point the overpressure had become 3 
percent higher than the strong shock prediction. The 
addition of 1 atmos to an inverse cube term gives a 
relation valid over a greater range of overpressures, 
deviating from the calculated curve by less than 5 
percent at 5 atmos. 


AP,=0.1567A,;3+1 atmos. (16) 


The shock radius (A,) is in the dimensionless units 
(energy/pressure-reduced) described in Sec. II. 
At still lower pressures the following empirical fit 
applies. 
0.137 0.119 0.269 


P= + —+ —0.019 atmos, 
rd dr? re 


for 0.1<AP,<10, 


or 0.26<,<2.8. (17) 


For all the calculations begun with a strong-shock 
the solid curve in Fig. 1 presents the peak overpressure 
as a function of shock radius. 
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Fic. 3. Peak or shock dynamic pressure (Q= 4pu*) versus shock 
radius for the point-source solution. The dynamic pressure is in 
atmos and the radius in units of (Etot/Po)t. 


The dashed curves represent the peak overpressures 
resulting from the hot isothermal spheres with normal 
density inside (beginning at 2002 and 122 atmos 
overpressure). Notice that the overpressure becomes 
indistinguishable from that of the strong shock at a 
radius where the mass of air engulfed by the shock is 
10 times the initial mass. In the strong-shock region 
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this corresponds to a drop in overpressure to about 17 
percent of the initial value, and in any case to a radius 
a little more than twice the initial radius. 

The dotted curve represents the overpressure from 
an initially cold isothermal sphere (normal tempera- 
ture, high pressure). Because of the slower rarefaction 
speed the shock pressure does not rise to the point- 
source value until quite late, i.e., until the rarefaction 
which moves inward from the initial pressure front has 
reached the center. 

In each case of an initial isothermal sphere there are 
some slight oscillations about the point-source over- 
pressure curve which are caused both by rarefactions 
and small shocks which form behind the front shock, 
move inward and reflect off the origin and then move 
out to overtake the shock front. (See a later discussion 
of these special features of the isothermal sphere 
problems.) 


Peak Dynamic Pressure (Q,) and Peak Particle 
Velocity (u,) versus Shock Radius (,) 


The particle velocity and the peak pressure at the 
shock front follow precisely the Hugoniot relation, 


SAP, 


gaa (18) 
(49+-42AP,)! 


WEL A A 
HIE 


4979 


6229 


.@) 500 R_ '!000 1500 
° 


Fic. 4, Pressure in atmos as a function of the Lagrange or mass position for the point-source solution at times indicated. The 
position is in arbitrary units of (Eto:/Po)?/1627.2, and the time is in units of (Etot/Po)#/Co when Cp is the pre-shock sonic velocity. 
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Fic. 5. Particle velocity in units of pre-shock sonic velocity as a function of Lagrange or mass position for the point-source solution 
at the times indicated. The position is in units of (Etot/Po)?/1627.2, and the time is in units of (Etot/Po)#/Co. 


The particle velocity follows its strong-shock de- 
pendence on shock radius down to as low an over- 
pressure as 75 atmos (Fig. 2); in fact an adequate fit 
over the entire range is given by the following: 


U,=0.30A 5-73. (19) 


This simple form is reflected in the peak dynamic 
pressure (0,=4p,0,") which falls somewhat faster than 
an inverse cube law down to equally low shock pressures 
(Fig. 3). The slow decrease in shock density (from 6 
for strong shocks to one for weak shocks) is responsible 
for the decay being steeper than the inverse cube of 
the shock radius. 


Pressure (P) versus Lagrangean Distance (Rj) 


The pressure behind the shock wave, shown in Fig. 
4 as a function of the Lagrangean (or initial position) 
variable, retains the strong-shock form until quite 
low pressures. Note, for example, that the ratio of 
central pressure to shock pressure remains 37 percent 
down to 20 atmos and decreases slowly below that to 
33 percent by 3 atmos. Beyond 3 atmos a negative 
phase develops with the pressure falling to as low as 
0.8 atmos near the center. 

In Fig. 4 the Lagrangean position is given in arbitrary 
units (Ro), related to the dimensionless unit (Aq) by a 


constant multiplier, 
Ro 


: 20 
1627.2 i 


To= EAg= € 


It may help in visualizing dimensions to assume Ry to 
be in centimeters; then the corresponding blast energy 
will be about that of 200 lb of TNT (4.210% ergs). 


Particle Velocity (uw) and Density (0) versus 
Lagrangean Distance (Ro) 


Figures 5 and 6 indicate the progression of particle 
velocity and density as functions of the Lagrangean 
position. Again the strong-shock form is dominant 
until as low as 3 atmos of shock overpressure. 

The particle velocity transforms gradually from its 
very linear form to one much like the overpressure at 
large distances as the shock wave goes from strong 
to weak. 

The density, which in the strong shock is zero at the 
origin, implying an infinite temperature, remains zero 
there since no mechanism is included (no conduction 
or radiation) for dissipating this high temperature. 

The density profiles in Fig. 6 represent strong-shock 
initial conditions at 200 atmos shock pressure. The 
dip in the curves that sit at the same mass point 
(Ro=150) is due to the sudden inclusion of a finite 
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Fic.'6. Density in units of pre-shock densityzas a function of Lagrange or mass position for the point-source solution 
at the times indicated. The position is in units of (Etot/Po)#/Co. 


atmosphere ahead of the shock, i.e., due to the abandon- 
ment of the strong-shock (infinite pressure ratio) 
assumption at that point. The Hugoniot relations give 
the shock density. in terms of the shock pressure as 


_GEDP+(y=1)_ 6PH1 
, (y—-1)P+ (y-+1) P+6 


For a strony shock (P=) the density ratio is six 
(o=6), but at two hundred atmos (P= 200) the ratio 
is only 5.83. 

The strong-shock solution is derived on the assump- 
tion that the atmosphere ahead of the shock wave 
has negligible effect. The effect, when non-negligible, 
is to raise the temperature through the shock to a 
higher value than that given by the strong-shock 
solution, i.e., a finite shock is hotter than would be 
predicted by the “‘strong-shock” theory. 


for y=1.4. (21) 


Gas Variables as a Function of Eulerian 
Position (2) 


The pressure, density, compression, and particle 
velocity are shown relative to their peak or shock values 
as functions of the Eulerian position (A) in Fig. 7. 
The strong-shock form dominates the first two sets, 
while the later ones show the characteristic positive 
phase followed by a longer, weaker negative phase 


and eventually by a return to near pre-shock values at 
the origin. 


Positive and Negative Phases versus Distance 


The durations of the positive phases for pressure and 
particle velocity are shown in Fig. 8. Although these 
durations should approach the same value at large 
distances, they still differ by 7 percent at a distance of 
=3.0(AP,=0.09). The figure indicates only the values 
for the point source solution. 

The duration of the negative pressure phase (D,-), 
unlike that of the positive phase, is nearly independent 
of distance, and has an average value of 1.22. The 
negative durations suffer somewhat from loss of 
numerical significance at late times in the calculation. 


Time Dependence of Pressure 


In Fig. 9, the curves of pressure versus time at various 
distances are given in units of the peak overpressure 
and the positive duration in order to illustrate the 
change in the rate of decay behind the shock. At 
increasing distances the drop in pressure is both 
relatively and absolutely slower since the positive 
phase duration is also increasing. Again, these curves 
are for the point-source solution. 

An exponential form approximates the time de- 
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Fic. 7. Hydrodynamic quantities in units of their peak values as a function of Eulerian or space position at the times indi- 
cated. The radius is in units of (Etot/Po)t. The solid curves represent overpressure (AP/AP,), the long dash curves represent 
particle velocity («/u,), the dot-dash curves represent density te/o), and the short dash curves represent compression [ (o—1)/ 
(ps—1)]. The shock values of the overpressure (AP,) for these times are 121.5, 20.10, 2.03, 1.01, 0.338, 0.0701 ordered on in- 
creasing time. For the particle velocity the peak values are 8.49, 3.37, 0.873, 0.524, 0.216, 0.0525, and for the density the peak 
values are 5.66, 4.68, 2.135, 1.625, 1.228, 1.0508. 


pendence of the pressure pulse: pressures less than one atmosphere (A>0.74), with the 
ep coefficient (a) specified by 
es (1—Z)e#4. (22) a=i+AP,, AP,<1. (23) 


For shock overpressures greater than one atmosphere 
Here Z is the time after shock arrival in units of positive the decay is not a simple exponential, since the early 
duration, AP, is the shock overpressure, and @ is portion requires a larger a than the later part. Allowing 
independent of Z. This form is satisfactory for over- the coefficient a to be a function of the time (2), it. 
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may be approximated by 
a=4+AP,[1.1— (0.13+0.2AP,)Z] 


for AP, less than 3 atmos, and for overpressures from 
3 up to 50 atmos by the form 


(24) 


b 
1+cZ 


a=a+ (25) 


where 
| —0,231+0.388AP,—0.0332AP2 for AP,<10 


0 for AP,>10 

AP,(0.88+0.072AP,) for AP,<10 
Sea for AP,>10 
c=8.714-0.1843AP,—104/(AP, +10). 


The time dependence of the pressure in the negative 
phase may be approximated by the form 


AP=14AP_w(1—w)e“ for 0.1<AP,<200 (26) 


5 


Tt, time 


° 
o i] 2 
X, distance 


Fic. 8. Duration of positive phase for pressure (D,*) and 
particle velocity (D,*) versus distance (Eulerian) where time is in 
units of (Etot/Po)?/Co and distance is in units of (Etot/Po)!. 


where AP_ is the peak negative overpressure (see Fig. 
10), and where w is the time measured from the end 
of the positive phase in units of the negative phase 
duration (D,-). A more accurate fit would allow the 
exponent to decrease to zero as the shock strength 
goes to zero. 


Dynamic Pressure versus Time 


Of interest also is the form of the dynamic pressure 
(Q=4pu*). Time plots of this function appear in Fig. 11. 
These curves are also normalized, but with the positive 
duration of the particle velocity (D+) (Fig. 8), and 
the peak value of the dynamic pressure (Q,) (Fig. 3). 

A similar exponential form approximates the dy- 
namic pressure, 


Q/Q.= (1—-Z)eP?, (27) 


where 8 is independent of the time (Z) for shocks of 
less than one atmos peak overpressure (AP, <1). The 
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Fic. 9. Pressure as a function of time where the pressure is in 
units of the peak pressure (AP,) (Fig. 1) and the time (Z,) is in 
units of the positive duration (D,*) (Fig. 8). The numbers indicate 
the corresponding peak overpressures (AP,). 


coefficient 8 changes with the shock strength, however: 

B=0.75+3.2AP., AP.<1. (28) 

Where the shock is stronger, a modification similar to 

that given in Eq. (25) for the overpressure is 
appropriate: 


p=d+f/(1+82), (29) 


50>AP.>1, 
{sett for AP,<3 


for 


—5.6+0.63AP, for 3<AP,<10 
0 for AP,>10 


d= 


f=6A40AP,, 
g=0.725AP,. 
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Fic. 10. Peak negative overpressure in atmos as a function of 
radial distance in units of (Ztot/Po)t. 
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Fic. 11. Dynamic pressure (Q=4pu*) as a function of time 
where the pressure is in units of the peak dynamic pressure (Q,) 
(Fig. 3), and the time is in units of the positive duration of the 
velocity (Du+) (Fig. 8). The numbers indicate the corresponding 
peak overpressures (AP,). This rapid decrease in effective pulse 
duration with increasing peak overpressure has the effect of 
making the higher dynamic impulses much the sharpest. 


These approximate forms agree with the numerical 
values of Q and AP in the positive phase to within 10 
percent (and over most of the range to less than 2 
percent) for values of AP, less than ten atmos. 


Positive Impulse 


The integrated positive overpressure (J,*) and the 
total positive drag pressure (/,,*) decrease with distance 
from the blast source in the manner shown in Fig. 12. 


Dot 
It= f AP()di, 
(30) 


Dut 
a 
0 : 


(Because of the units employed in this paper, [,* 
does not become the usual dynamic impulse until it is 
multiplied by y= poC0?/Po.) 

Both of these impulses may be fitted by simple 
powers of the radial distance for shock overpressures 
less than two atmospheres: 


Tt =0.043-, 
T.+=0.004).7, 


AP, <2 (31) 


Negative Impulse 


The negative overpressure impulse in the range 
below 20 atmos of peak overpressure can be expressed 
to within 5 percent by 


1,-=4A4P_D,-. (32) 
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Since the duration of the negative pressure is con- 
stant at Dp-—~1.22, 


Ip-0.61AP_. (33) 


For distances greater that \=0.5 the negative peak 
overpressure approaches zero inversely as the distance 


AP_~—0.086\, »>0.5, (34) 


and the negative pressure impulse goes approximately 
like 
Ip-0.052\,, A> 0.5. (35) 
Special Features of the Isothermal Sphere 
Problems 


In the gas dynamics resulting from the release of an 
initially static high-pressure sphere, an inward moving 
shock forms behind the rarefaction wave that first 
runs in from the surface of the sphere. This inward- 
directed shock was predicted by Wecken’ and discussed 
by McFadden,* Shardin® and others.4 It does not 
acquire a net inward velocity until the rarefaction 
has reached the center, but after that it moves in and 
reflects at the origin, and then races outward to eventu- 
ally overtake the main shock. 

This second shock grows from zero strength to 
presumably an infinite pressure ratio at the origin. 
On reflection, this shock moves outward, decaying in 
strength about as the inverse first power of its distance 
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Fic. 12. The radial dependence of the positive pressure impulse 
(7,*) the negative pressure impulse (Jp-), and the positive 
dynamic impulse J,* in units of atmospheres [where time is in 
the dimensionless units (Eto:/Po)t/Co] the radial distance is in 
units of (Ftot/Po)t. 


7 F, Wecken, Z. angew. Math. u. Mech. 30, 270 (1950). 

8 J. A. McFadden, J. Appl. Phys. 23, 1269 (1952). 

®H. Schardin, ‘Measurement of spherical shock waves,” 
Commun. Pure Appl. Math. (Institute of Mathematical Sciences, 
New York University) VII, 223 (1954). 


SPHERICAL BLAST WAVES 


from the origin. The general expansion and outward 
motion of the gas behind the main shock is responsible 
for this apparently modest decay rate. 

When this second shock encounters the outer surface 
of the gas that was initially inside the sphere, a trans- 
mitted shock continues out and a reflected shock is 
sent inward. The transmitted shock overtakes the main 
shock and increases it by as much as 20 percent. 
Progressively weaker shocks follow from the reflections 
at the interface and again at the origin. This repeated 
shocking of the gas near the origin changes the tem- 
perature profile at later times from its initial isothermal 
nature to nearer the point-source distribution which is 
characterized by a high central temperature falling off 
rapidly with radius. 

Some detail is lost in the vicinity of the origin since 
the shocks are always spread over a number of mesh 
points. In one case the problem was rerun with mesh 
sizes about one-fourth the original size, and the inward- 
moving shock showed some appreciable discrepancy 
near the origin. However, this discrepancy may be 
attributed mainly to the difficulty in identifying the 
shock when the rounding is comparable to the shock 
radius. After reflection the difference vanishes again 
in spite of the quite different histories near the center. 


IV. CONCLUSIONS 
Point-Source Solution 


At increasing distances from a finite but sudden 
source of energy, the resulting blast wave will appear 
more and more like that from a point source. The 
blast resulting from an initial isothermal sphere of gas 
at rest will assume the general shape and values of the 
point-source solution (to within 10 percent) after the 
shock wave has engulfed a mass of air 10 times the 
initial mass of the sphere. Prior to this, the shock 
strength is less than that of the point-source shock, and 
the inward-traveling rarefaction has not reached the 
center. 

A point source should leave a higher temperature 
and consequently a larger percentage of energy near the 
origin. This energy, no longer available to the shock 
wave, should effect a reduction in the shock radius 
(for a given overpressure). However, no appreciable 
difference in the low-end of the shock overpressure- 
radius relation appears on comparing the point source 
and isothermal sphere solutions. In fact, perhaps 
because of multiple shocking of the inner regions in the 
isothermal sphere problems, the distribution of residual 
energies (per unit volume) and pressures are nearly 
identical around the origin at a time when the shock 
has progressed to 6 or 7 times the initial radius. 
Although temperatures will remain different, since the 
point source has infinite temperature at the origin 
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{and zero density), the average temperature (or 
density) of that gas initially inside the isothermal 
sphere approaches (within 10 percent) the average 
temperature for a corresponding mass around the point 
source, 

For a source of considerable initial mass, the peak 
pressure may become quite small before the shock 
has engulfed a mass of gas 10 times larger, and the 
blast wave may, therefore, remain quite different from 
the point-source solution throughout regions of interest. 
Such is the case to some extent with high explosives 
where initial charge shapes will influence the blast wave 
at all significant pressures. It is true to an even greater 
extent in the spherical equivalent of a shock-tube type 
of blast, where a gas at high pressure but at normal 
temperature is suddenly released, as in the normal 
temperature isobaric sphere problem (dotted curve 
Fig. 1) described in Sec. I. 


Application to Blast in Air 


The ideal gas assumption is reasonably valid in air 
for shock pressures less than 10 atmos. Above that the 
gamma in the expression for the internal energy 
[Eq. (7) ] ranges down to as low as 1.133 and up to the 
monatomic value of 1.667. But since the mass of gas 
which has experienced shocks stronger than 10 atmos 
is a small part (5 percent) of the engulfed mass by the 
time the shock overpressure is down to 1 atmos, the 
blast wave for constant gamma should be reasonably 
correct for air at most interesting shock pressures. 

In the case of shock wave problems involving only 
one space variable, the artificial viscosity technique 
for numerical integration appears to be very satisfactory 
over large ranges of pressure and entropy change. 
Unfortunately, it cannot be expected to yield details of 
a shock wave impinging on a singular point such as the 
origin in spherical or cylindrical geometry, since the 
nature of this method is such as to spread the shock 
front over a number of grid points. 
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